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We apply the background field (BF) method to Non-Relativistic QCD (NRQCD) on the lattice 
in order to determine the one-loop radiative corrections to the coefficients of the NRQCD action 
in a manifestly gauge-covariant manner by matching the NRQCD prediction for particular on-shell 
processes with those of relativistic continuum QCD. We explain how the BF method is implemented 
in automated perturbation theory and discuss the technique for matching the relativistic and non- 
relativistic theories. We compute the one-loop radiative corrections to the a ■ B and Darwin terms 
for the NRQCD action currently used in simulations, as well as the one-loop coefficients of the spin- 
dependent 0(a 2 ) four-fermion contact terms. The effect of the corrections on the hyperfine splitting 
of bottomonium is estimated using earlier simulation results [l]; the corrected lattice prediction is 
found to be in agreement with experiment. Agreement of the hyperfine splitting of bottomonium 
and the _B-meson system is confirmed by recent simulation studies [2J [3] which include our NRQCD 
radiative corrections for the first time. 
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INTRODUCTION 



The decays of hadrons containing bottom quarks pro- 
vide some of the most stringent tests of the Standard 
Model (SM) so far, but direct simulations of the bottom 
quark in lattice QCD remain fraught with problems due 
to discretization effects arising from the large mass of the 
bottom quark in conjunction with current limitations on 
achievable lattice spacings. An alternative is provided 
by Non-Relativistic QCD (NRQCD) [4], an effective field 
theory for heavy quarks whose use in describing heavy- 
flavour hadrons has so far met with considerable success 
PQ and which represents one important approach to ob- 
taining accurate predictions for flavour physics observ- 
ables and testing the limits of the SM. However, until 
recently the NRQCD actions used in lattice simulations 
did not include radiative corrections to the action, lim- 
iting the accuracy to which such important quantities as 
the T-r/b hyperfine splitting could be predicted [5] . This 
is in stark contrast to the situation of Non-Relativistic 
QED (NRQED), for which the radiative improvements to 
the action have long been known, leading to highly pre- 
cise theoretical predictions for muonium hyperfine struc- 
ture and for positronium decay [SHE]- Achieving sim- 
ilar precision for NRQCD predictions by including the 
radiative improvements in the NRQCD action is there- 
fore highly desirable. However, a crucial difference to the 
NRQED case is that the strongly interacting non-abelian 
nature of QCD and NRQCD imposes confinement and 
calls for a lattice implementation of NRQCD; this means 
that the complete l/(Mo)™ structure of all quantities 
must be retained, as opposed to the situation in (con- 



tinuum) NRQED, where there are ways to drop terms 
in (A/M) n consistently. A further complication arises 
from infrared (IR) divergences, which turn out to play a 
significant role in the non-abelian case. 

In this paper, we follow up on our letter (9] where 
we presented the first calculation of radiative correc- 
tions to the lattice NRQCD action using the background 
field (BF) method. We proceed by computing the effec- 
tive action in both continuum QCD and lattice NRQCD 
at one loop and matching the latter term by term to the 
non-relativistic reduction of the former. In particular, 
we determine the 0(a s ) correction to the coefficient C4 
of the chromomagnetic cr ■ B operator and the leading 
contributions to the coefficients of the spin-dependent 
four-fermion contact operators in the NRQCD action, 
as required in order to enable a more precise determi- 
nation of the hyperfine structure of heavy quarkonia in 
NRQCD. Using our results, we are able to estimate the 
0(a 2 s ) correction to the hyperfine splitting of the S-wave 
bottomonium states using simulation data from [1], giv- 
ing a corrected value of 72(3) (5) (3) MeV, which agrees 
with the experimental value of 69.3(2.8) MeV [TP] .The 
correction has the additional beneficial effect of reduc- 
ing the lattice spacing dependence, placing the remain- 
ing 0(a 2 ) discretization errors well below other sources 
of error. 

We remark that our calculation reported here com- 
pletes the one-loop radiative improvement of quark bilin- 
ear terms to 0(v A ) in the NRQCD action for which the 
improvement of the purely kinetic terms was reported in 



0: see section III 



In section [H] we discuss in some more detail the ap- 
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proach to matching using the BF method and discuss 
the gauge invariance of our results; in sections IV andfV} 
we proceed to apply the BF method to the matching of 
the chromomagnetic and Darwin terms in the NRQCD 
action. The four-fermion interactions in NRQCD are the 
topic of section |VI[ where we derive results for the spin- 
dependent four-fermion operators and explain the diffi- 
culties in matching their spin-independent counterparts, 
which will be the topic of a future paper. We apply our 
results to the hyperfine splitting and S-wave mass shifts 
in section |VII| before summarizing our results in section 

ehd 

In the following we denote the perturbative expansion 
for a generic quantity w as w — ^2 n=0 w^a™. Where 
contains IR divergences, the IR finite part is denoted 



w 
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with a tilde: 



II. THE BACKGROUND FIELD METHOD FOR 
LATTICE NRQCD 

The BF method [TTHT5] is well established as a tool 
to compute the effective action in quantum field theory, 
which has a number of very attractive features: firstly, 
QCD in background field gauge (BFG) satisfies a set of 
abelian-like Ward identities reducing the number of cal- 
culations needed in order to renormalize it. Secondly, the 
existence of a residual gauge invariance in BFG implies 
that only gauge-covariant operators can appear in the 
effective action, which is of particular importance when 
considering operators of dimension D > 4 as required 
in an effective theory, where a loss of gauge-covariancc 
would herald a proliferation of additional operators. 

A method to derive a unique effective action is given by 
the Vilkovisky-DeWitt (VDW) technique and one choice 
is to use the Landau-DeWitt gauge corresponding to the 
BF generalization of Landau gauge [TB] . Whether or not 
this approach is applicable to our matching procedure 
remains to be investigated. However, for our purposes, 
the VDW method is not necessary since we are able to 
perform our matching solely using on-shell quantities in- 
cluding S-matrix elements which can be given an unam- 
biguous physical interpretation. 

The BF method is indispensable to our matching pro- 
cedure. As an effective theory, NRQCD has operators of 
dimension D > 4 appearing in the action, and BFG is 
needed to ensure a gauge-covariant form of the effective 
action; this is not guaranteed nor likely without using 
the BF formalism. It should be noted that the appear- 
ance of gauge-non-covariant D > 4 operators in the effec- 
tive action would not in and of itself be incorrect; how- 
ever, by obscuring the gauge symmetry of the theory, 
they would hinder a physical interpretation of the results 
and significantly complicate the calculations. Moreover, 
without the use of BFG, a further source of complexity 
would arise from the appearance of ultraviolet (UV) log- 
arithms in the coefficients of operators in the effective 
action, whose contributions to physical processes would 



have to cancel against the contributions from the addi- 
tional non-gauge-covariant operators. 

Finally, the BF method makes our work easier and 
tractable because the QED-like Ward identities in BFG 
are enough to render the one-particle irreducible (1PI) 
vertex function finite, so that only the gluonic self-energy 
renormalizes the coupling, whereas the BF itself is not 
renormalized. Since these statements hold in both QCD 
and NRQCD, we are able to match the two theories by 
equating two quantities that are UV finite and thus we 
can use different regulators on either side. In particu- 
lar, the QCD vertex can be calculated analytically in 
the continuum using dimensional regularization (DR), 
or numerically on fine lattices approaching the contin- 
uum limit, where the latter is particularly convenient for 
checking the gauge-parameter independence of the result 
since analytical calculations for arbitrary values of the 
gauge parameter tend to become rather involved. 

The remaining issue is how to regularize the infrared 
(IR) divergences which arise at intermediate steps in our 
calculation. The radiative corrections to the coefficients 
in the NRQCD action are, and must be, IR finite: they 
cannot depend on any scheme to regularize IR diver- 
gences. Two kinds of IR divergences arise. The first 
are those that occur in the continuum calculation of the 
diagrams. These divergences are a-priori independent of 
the non-zero lattice spacing a and must match directly 
between the relativistic QCD and NRQCD calculations. 
The second are lattice artifact IR divergences that occur 
solely in the NRQCD one-loop calculations, and the con- 
tributions from the set of diagrams under consideration 
must cancel consistently within the NRQCD calculation 
since they depend on a. That the IR divergences match 
or cancel in this way is a strong consistency check on our 
calculations. In NRQCD we are also able to calculate 
both kinds of IR divergence analytically in every case 
considered and also numerically and agreement between 
the two approaches is a further check on our results. 

Whilst it might, in principle, be possible to eliminate 
the continuum IR divergences by subtracting the inte- 
grand in continuum QCD from the equivalent one in 
NRQCD, this approach is numerically very difficult. It 
is compounded by the need to also eliminate the lattice 
artifact IR divergences in the same way, otherwise the 
use of an IR regulation scheme is, in any case, inevitable; 
this would be numerically very hard to accomplish. We 
cannot use dimensional regularization to regulate the IR 
divergences since these are not applicable for lattice field 
theories. There are a number of other methods that can 
be used in our case. 

Twisted boundary conditions (TWB) can be applied 
corresponding to matching the non-abelian gauge theo- 
ries in a region of finite extent L in two or three spatial 
dimensions [TTHT^] . The TWB eliminate the zero modes 
for the non-abelian gauge field and give a gauge invari- 
ant regulation of IR divergences with a mass-scale 2ir / 3L 
for the SU(3) gauge theory. Although there are extra 
complications in the background field approach, TWB 
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are relatively straightforward to implement for numeri- 
cal lattice calculations and have been used to carry out 
improvement calculations 19. 20 . However, they are dif- 
ficult to implement analytically for the continuum QCD 
calculations and although there is progress in this area 
we choose not to implement the TWB. 

The NRQCD action is a derivative expansion and the 
strategy is to match the expansion in independent ex- 
ternal momenta for each amplitude considered. IR di- 
vergences can be regulated by using a small but non-zero 
momentum transfer, q say, in the diagram concerned. Af- 
ter cancellation of the IR singularities we may safely set 
q = 0. The application of this method to diagrams with 
more than one loop needs careful thought but is certainly 
an option for one-loop calculations. 

We regulate the IR divergences in both QCD and 
NRQCD using a gluon mass fi. This is known to be 
correct for diagrams which have a QED-like topology: in 
QCD the difference is a simple overall colour factor. In 
general, a non-zero gluon mass in QCD breaks gauge in- 
variance and in our calculation the one-loop diagrams for 
matching the quark bi-linear operators, which give the 
chromodynamic form factors, contain three-gluon ver- 
tices and the use of a gluon mass IR regulator needs 
justification. We note that at one-loop the diagrams are 
identical to those that would arise in the Curci-Ferrari 
theory [5TJ [55] which is renormalizable and known to re- 
cover QCD in the [i — > for gauge-invariant quantities; 
the introduction of a non-zero gluon mass by extending 
QCD to the Curci-Ferrari formulation is a mechanism 
for regulating IR divergences [2TH2T)] . We match on-shell 
gauge-invariant physical processes and in our final result 
all IR divergences cancel to give IR-finite gauge-covariant 
counterterms in NRQCD and the limit /i — » can be 
taken. Diagrams with more than one loop in Curci- 
Ferrari theory will contain extra contributions compared 
with the same process in QCD. The Curci-Ferrari theory 
on a lattice is formulated and discussed by von Smekal 
et al. [30]. 

For the case of matching the quark bi-linear oper- 
ators another approach is to introduce gluon masses 



by the spontaneous symmetry breaking of SU(3) — > 
(7(1)3 x U(l)s [31] and calculate in the renormalizable 
i?£ gauge. The unbroken abelian gauge groups are gen- 
erated, respectively, by the T 3 and T 8 SU(3) genera- 
tors. In general, QCD is not recovered in the zero mass 
limit since there remain massless scalar fields. However, 
the one-loop diagrams that we calculate are identical in 
both QCD and the spontaneously broken theory. In this 
case, our present calculation corresponds to computing 
the quark bi-linear form factors for the U(l)s (massless) 
gluon. 

An alternative justification was established by Kaut- 
sky [32]. In BFG the Ward Identities in QCD give a 
QED-like relationship between the quark wavefunction 
renormalization constant, Z2, and the vertex renormal- 
ization constant, Z\, namely Z\ — Zi\ a violation of this 
equality would signal a breakdown of gauge invariance. 
Kautsky showed in BFG that at one loop this equality 
does hold when the IR divergences in the diagrams con- 
cerned are regulated by a non-zero gluon mass. This is 
a crucial test of this approach for regulating IR diver- 
gences. 

An important observation, that we verify, is that the 
results of calculating the physical process used for the 
matching procedure in NRQCD and QCD are separately 
independent of the gauge parameter multiplying the BF 
gauge-fixing term in the action: the NRQCD radiative 
counterterms are both gauge-covariant and independent 
of the gauge parameter. 



III. THE NRQCD ACTION 

The lattice NRQCD action up to and including 0(v 6 ) 
operators is given by [4] 

Snrqcd = ^2^(x,t) [ip(x, t) - K(t)iP(x, t - a)] , 

(1) 

with the kernel 



The quark Green function satisfies the evolution equa- 
tion 

G(x, t + a) = K(t + a)G(x, i) + S Xt0 5 t+afi , (3) 

where G(x, t) vanishes for t < 0. The kernel K is con- 
structed so as to give an evolution that is symmetric with 
respect to time reversal, and leads to a smaller wave func- 
tion renormalization than some other formulations [33] . 
The parameter n was introduced to prevent instabilities 



at large momenta due to the kinetic energy operator |34j 
fit is easy to see that the Green function defined by Eq. 
M diverges if aH l ™ pr /2n > 2). For a given /3 and lattice 
spacing we therefore have a minimum n value. For the 
values of lattice spacing used by the HPQCD collabora- 
tion n = 2 or n — 4 suffices. 

The form of the evolution for the heavy quark Green 
function with a single time derivative leads to an initial 
boundary value problem that is relatively easy to solve, 
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whereas introducing further terms containing higher- 
order temporal derivatives would create a computation- 
ally much more demanding boundary value problem. We 
therefore do not improve the discrete temporal derivative 
in the usual way but instead modify the kinetic energy 
operator by using 
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(4) 



The unimproved forward and backward derivatives are 
defined by 

(A+f)(x) = ±(U fi (x)f(x + fi)-f(x)) ) 

(A- /)(.*) = !(/(*) - U^{x)f{x - A)) , (9) 

and the unimproved higher order lattice derivatives are 

3 



where Hq is the standard NRQCD kinetic operator, Hq = 
A&/2M. 

To understand this correction consider the Green func- 
tion of the standard NRQCD kinetic operator in the ab- 
sence of gauge fields, 



A (2n) = ^ (A (+) A (-) ) , 
3=1 



(10) 



The 0(a p )-improved symmetric derivative is given by 



G(x,t + a) = 1- 



aH, 



2)1 



A(±) 



A 



2n 



G(x,t) 



(±) _ ^a (+) a (±) aH 

g J 3 3 



(11) 



(5) 



which for a single time step has an exact exponential form 
in terms of 



B and E are the improved chromomagnetic and chro- 
moelectric fields, constructed from standard cloverleaf 
operators, and improved to 0(a 2 ): 



2n aH 



r> _ __, pimp 
J^k — r, t klm- r i m , 



E, 



rpimp 
^k4 ' 



(6) 



pimp _ ^ rp 
1 p,v ~ o f 



so that we can remove the leading 0(a) discretization 
error by introducing H % ™ pr . 

The interaction terms for the "full NRQCD" action are 
given by, 



Fp,i> — „ (Ifiv i)i 
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(12) 
(13) 

I, 

(14) 
(15) 
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A(±) • E - E • A (±: 



J2 U a (x)U p {x + a)U- a (x + a + f3)U-p(x + p) , 

(16) 



c 4 — cr • B - c 3 -|W°- ' I A (±) x E - E x A^ 

2M 8M 2 
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2AM 



(A^ 2 )) 2 
16nM 2 



C6 



(7) 



including operators with coefficients C5 and cq to re- 
move the leading discretization artifacts in the improved 
NRQCD kinetic operator at O(o 4 /). The operators have 
been normalized such that a — 1 at tree level. The one- 
loop radiative corrections to c\ , C5 and cq have been re- 
ported in [2]. We can also add further spin-dependent 
interaction terms to obtain the "full spin v 6 NRQCD" 
action given by 



5H, 



-fi 
-h 



3l S f A(2 ) _ ^xE-ExA^ 



64Af 4 
2 



{a( 2 \<x. (z 
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+ f 3 ^(T • E x E . 

Ji 8M 3 



)} 

(8) 



where again the operators have been normalized such 
that at tree- level fi = 1. 



where {(a,/?)}^ = (//, v), {v, -fi), (-fj,, -v), (-v,fi) for 
A* 7^ 

Mean-field (tadpole) improvement is implicitly in- 
cluded in the definition of the links by replacing {/„ t— > 
Ufi/(uo), where uq is the Landau gauge mean link. Since 
the gauge links are unitary, some of the correction terms 
in F lmp will be only four, rather than six, links long due 
to cancellations of the form U^U^, and so should carry 
a factor of \ (rather than \). This is implemented by 

u Q u 

the addition of a further correction term |35j and must be 
taken into account when calculating mean-field improve- 
ment contributions. 



IV. CHROMOMAGNETIC INTERACTION 

The hyperfine splitting of S-wave bottomonium, Mx — 
M Vb , provides a high-precision test of NRQCD [5 . The 
size of the hyperfine splitting is expected to be ap- 
proximately proportional to the square of the coefficient 
c 4 of the chromomagnetic operator in the NRQCD ac- 
tion. While most NRQCD simulations so far have set 
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(a) 



(1>) 






(e) 



FIG. 1. Vertex correction Feynman diagrams. Referred to as 
(a) QED (b) nonabelian (c)&(d) swordfish (e) algae (f) ankh 
diagrams. 



C4 = 1, this coefficient will receive radiative corrections 
at 0{a s ), which have the potential of affecting the pre- 
dicted value for the hyperfine splitting significantly, so 
that without a determination of these radiative correc- 
tions, large systematic errors have to be included in any 
NRQCD calculation of the hyperfine splitting. For ex- 
ample, in [1] the ground state hyperfine splitting is pre- 
dicted to be M r - M m = 61(U)MeV, where the domi- 
nant errors are from the missing 0(a s ) corrections to C4. 
This prediction is smaller than the experimental value of 
M r -M Vb = 71.4(+2.3)(-3.1)(3.7), [35] although within 
the large theoretical errors, theory and experiment agree. 



At the one-loop level, the contributions to the effective 
action relevant for determining the chromomagnetic op- 
erator arise from the vertex diagrams shown in Figure [T] 
As the QCD analogue of the magnetic moment operator 
in QED, the chromomagnetic moment operator consti- 
tutes the leading contribution to the hyperfine splitting 
within hadronic states, and by tuning of the coefficient C4 
of the chromomagnetic operator in the NRQCD action, 
we can match NRQCD to QCD as far as the hyperfine 
splitting is concerned. 



In this section we compute the radiative improvement 
correction to C4 which takes the form C4 = l + a s c^ with 
C4 = A+B log(Ma), where a is the lattice spacing, M is 
the heavy quark mass, and A and B are constants which 
we calculate. In continuum QCD, all corrections can be 
computed analytically following standard techniques; for 
the calculation in lattice NRQCD, we employ the auto- 
mated Feynman rule packages HPSRC and HiPPy |37| . 



A. Continuum QCD Calculation 

The QCD effective action contains the following terms 
bi-linear in the fermionic fields, 



- G^ v F 

2m 



, (17) 



where, because of BFG invariance, the matter terms and 
\t A^> vertex are packaged into a gauge covariant deriva- 
tive so as to obey the Slavnov-Taylor identity Z\p = Z^ 1 ■ 
The constants Z are thus related to the usual chromo- 
dynamic form factors via Z\p = -Fi(O) and 5Z a — ^2(0). 
The first term is multiplicatively renormalized to 



R- 



(18) 



and BFG ensures that A is not renormalized. The second 
term rcnormalizcs to 



2m R 



with 



Z„ Zo z n 



(19) 



(20) 



Since QCD is renormalizable, the absence of a Pauli term 
in the bare QCD action ensures that b a is UV finite; it 
moreover implies that Z a = 0(a s ), so that to leading 
order we can set Z 2 = Z m = 1. 

After performing a non-relativistic reduction using a 
Foldy-Wouthuysen-Tani (FWT) transformation [381139] . 
Eqs. (fl8|) and (Il9| reduce to 



B 



t a 

R 2m R 



(21) 



where B is the chromomagnetic field. It is worth not- 
ing the relationship between various notations for these 
renormalization constants: at the one loop level 

1 + b a = 1 + Z a Z 2 = (Fi(0) + F 2 (0))Z 2 - G M (0)Z 2 . 

Note that there is no factor of Z m in these expressions, 
since the Gordon decomposition of the tree-level ^A^ 
term is between on-shell spinors and so we automatically 
have the renormalized mass in the resulting chromomag- 
netic term. As stated before, b a is UV finite; this is cru- 
cial for our analysis, since it enables us to directly equate 
results obtained on the lattice to those obtained in the 
continuum since the difference between the schemes for 
UV regulation is then irrelevant. 

The two diagrams contributing to & CT are given in Fig- 
ure [l] (a) and (b), and we introduce a small gluon mass 
/x as an IR regulator. A straightforward analytic calcu- 
lation then gives 



3a s pi 13a s 
ba = -7T- log — + -r— 

ZlT 771 D7T 



(23) 



at one loop. This continuum calculation was verified us- 
ing the HiPPy and HPSRC packages by reproducing the 
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same result numerically in the limit ma — ¥ 0; in this way 
it was also confirmed that the result was gauge-parameter 
independent. Note that the choice of BFG was crucial to 
ensure a UV finite result. 



B. NRQCD Calculation 

The leading spin-dependent term in the effective action 
for NRQCD is 



r NRQCD Z NR^t^l% + 



(24) 



which, after renormalization, becomes 



pNRQCD = CiZTzrzT ^°_B^ R (25) 

ZM R 

Since the chromomagnetic operator is present in the 
NRQCD action at tree level, Z a is of the form = 
1 + Z„ R '^a s . Requiring that the anomalous chromo- 
magnetic moment to be equal in QCD and NRQCD, we 
obtain the matching condition 



„ 7NR7NR7NR 



i + z a , 



which at tree level and one-loop order yields 



,(o) 
-4 

„(i) 



= 1, 



(26) 



(27) 



^NR,{1) _ ^NR,(1) _ ^NR,(1) 



Note that mass renormalization has to be included in 
NRQCD since the chromomagnetic operator is now in- 
cluded at tree-level. 

All of the renormalization constants involved are UV fi- 
nite and consist of two distinct contributions: besides the 
ordinary diagrammatic contributions (which we denote 
simply as Z), there are also contributions from mean-field 
improvement (Z tad ). The NRQCD dia grams contribut- 
ing to Za R ^ are shown in Figure [l] (a)-(e). Note that 
diagrams (c)-(e) receive contributions not only from lat- 
tice artifacts, but also from higher-adicity vertices that 
are also present in continuum NRQCD. 

Having once coded the Feynman diagrams in Figure [I] 
using the HPSRC package, we can repeat the numerical 
evaluation of these diagrams for each NRQCD action of 
interest by using the HiPPy package to automatically 
generate Feynman rules for that action along with the 
Symanzik-improved gluonic action. 

By using the residual gauge-invariance of the effective 
potential in BFG, the contribution to the negative of the 
effective action from the diagrams can be restricted to 
the form 



-z^ (1) ^(P+ q )' iei3k ^ Mq) m 



, (28) 



where i — 1,2,3 and terms of higher order in the exter- 
nal momentum q have been omitted. Without BFG, non- 
gauge invariant operators would be generated as counter- 
terms, greatly increasing the difficulty of the calculation 
and the subsequent use of the improved action in simu- 
lations. Using the automatic differentiation [4012T] and 
spinor manipulation facilities built into the HPSRC pack- 
age, the coefficient of interest can be isolated as 



<:) 



Mig-Tr (5T^ A2 (p,0)a 3 ) = gZ* 



(29) 



where p = 0, with po taken such that the quark is on- 
shell using the lattice equation of motion as implemented 
in the HPSRC package. 

We integrate over the temporal component ki of the 
loop momentum k by contour integration over the unit 
circle w — e tki in the complex plane. The positions of 
the poles in the integrand arising from the gluon and 
quark propagators in the Feynman diagram are functions 
of the spatial loop momentum k. To ensure that no poles 
cross the contour of integration as k varies, we change the 
contour to be a circle of radius r, choosing r such that 
the contour lies exactly half-way between the outermost 
interior pole and the innermost exterior pole. Since there 
are no physical intermediate states and thus no branch 
cut in any of the diagrams this strategy is always possible 
|42) . Since it is known that the poles of full NRQCD and 
improved gluonic actions always lie further away from 
the contour than their unimproved counterparts [43] . it 
is therefore safe to reduce the effort by locating the poles 
corresponding to the unimproved action and shift the 
contour appropriately. 

The renormalization constants 
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1 + Re£(0) +Im 
9 2 E 



Z* R,(1) = 1 + MRe 



dp 2 . 



dp 4 
- Im 



p=0 

dpi 



(30) 
(31) 



are determined from the quark self-energy 
—Y^{—p,p) — £(p), which is given by the diagrams 
shown in Figure [2j 

Both Z^'^ and Z& contain logarithmic IR diver- 
gences. For our purposes, it sufficient to evaluate their 
sum for which the IR logarithm is known from analytical 
results. We can therefore determine the IR finite dia- 
grammatic contribution Zx from fitting our numerical 
results with the form 



Z NR,(1) +Z NR,(1) = £NR,(1) +Z2 



NR,(1) 



2tt 



log /ia. (32) 



We evaluate the diagrams for 10~ 8 < fi < 10~ 6 such 
that any (na) 2 lattice artifacts can be neglected in the 
fit. The logarithmic IR divergence combines with the IR 
logarithm from the QCD result above to yield an overall 
logarithmic contribution — log(Ma) to C4 . 
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(a) 



(b) 




FIG. 2. Quark two-point function Feynman diagrams. Re- 
ferred to as (a) rainbow and (b) tadpole diagrams. 



Z a and Z m also receive corrections from tad- 
pole improvement, whereas there are no tadpole correc- 
tions to Z™'^. Using a Mathematica [H] notebook 
developed in the course of earlier related work [42j 05] , 
the tadpole corrections are determined by symbolically 
substituting mean field corrections, U —> U /uq with 

(2) 

uq = 1 — a s u y , into the NRQCD action. The tadpole 
corrections are dependent on the details of the NRQCD 
action used, in particular also on the value of the stability 
parameter n. 

For the full v 4 NRQCD action we find 

~tad,(i) _ _ A , 3 ^ U W (33! 
Zm ~ U (Ma) 2 J ' 

The tadpole contribution to Z& arises from two 
sources: the application of mean-field improvement 
to the improved field-strength tensor, and the cross- 
multiplication of the tree-level a-B term with the tadpole 
corrections terms in Hq pQ . The overall result for the full 
v 4 NRQCD action is 

~tad,(l) = A3 J3_ _ _J5 3_\ (2) 

\3 AMa 8n(Ma) 2 4(Ma) 3 J ' 

We take the one-loop contribution to the Landau mean 
link to be u { 2) = 0.750 [46]. 



C. Results 

The final result for the radiative correction to C4 is 
given by 

„(1) _ 13 _ 7NR,(1) _ ~NR,(1) _ ~NR,(1) 
4 — 67r a 2 

- Zl ad ' (1) - Z* alL(1) - A log Ma. (35) 

Results for several different NRQCD actions, all used 
with the tree-level Symanzik gluon action, are summa- 
rized in Tables UllIIII The values of Ma were chosen to 



correspond to the lattice spacings for lattices used by the 
HPQCD collaboration pQ. 

We note that the tadpole corrections to Z a and Z m 
in each case very nearly cancel the diagrammatic contri- 
bution, demonstrating that they are working as intended 
by reducing the coefficient of a s in the perturbative se- 
ries. Without such a tadpole correction we find C4 ~ 4, 
meaning that the viability of the perturbative expansion 
would be highly questionable. 

On the other hand, the radiative corrections for the 
full v 4 NRQCD action with stability parameter for n — 2 
or n = 4 differ very little; in the case of the full spin 
v 6 NRQCD action it appears that for small values of 
Ma the correction increases slightly, but for larger values 
of Ma the corrections are very similar to the v 4 case. 
These results suggest that the radiative corrections for 
the chromomagnetic operator are relatively independent 
of the details of the NRQCD action. 

Note that for all actions and ranges of Ma the total cor- 
rection is positive: the constant part of the correction is 
larger than the (negative) logarithmic contribution. This 
refutes the claims of Penin [57], who does not include a 
calculation of the Ma-independent constant contribution 
to C4 which is responsible for being positive. 

In Figure [3] we plot the dependence of the radiative 
correction C4 against Ma for full NRQCD n = 4. While 
the expected divergence in the Ma —> limit can be 
seen, for values of 2 < Ma < 4 the correction varies only 
slowly, demonstrating that for the range of lattice spac- 
ings used by the HPQCD collaboration the perturbative 
improvement corrections are under control. 

1.8 . . . . . . . 1 

• 

1.6 
1.4 
1.2 

J? 10 ' 

"0*0.8 . * ' 

0.4 
0.2 

°8.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 

aM b 

FIG. 3. Lattice spacing dependence for Chromomagnetic op- 
erator correction, full 0(v 4 ) NRQCD n = 4 action. 



V. THE DARWIN OPERATOR 

The vertex correction diagrams of Figure [l] also renor- 
malize the D • E operator (conventionally called the Dar- 
win term in analogy with atomic physics) in the effective 



Ma 


1.95 


2.8 


4.0 


zP + zP 


-5.164(7) 


-4.913(6) 


-4.739(6) 


7(1) 


1.512(1) 


1.022(3) 


0.723(2) 


^tad,(l) 


4.387 


4.077 


3.841 


„tad,(l) 


-1.092 


-0.787 


-0.641 


c (1) 


0.728(7) 


0.799(7) 


0.842(6) 



TABLE I. Renormalization parameters for Chromomagnetic 
term, full NRQCD n = 2 



Ma 


1.9 


2.65 


3.4 


z<P+z<» 


-5.17(1) 


-4.94(1) 


-4.80(2) 


7(1) 


1.56(1) 


1.08(1) 


0.84(1) 


™tad,(l) 


4.43 


4.129 


3.95 


r7tad.fl) 


-1.12 


-0.82 


-0.69 


c (1) 


0.68(1) 


0.78(2) 


0.82(3) 



TABLE II. Renormalization parameters for Chromomagnetic 
term, full NRQCD n = 4 



action. The primary effect of this operator is to change 
the effective potential for the wavefunction at the origin; 
since only states with L = have a non-vanishing wave- 
function at the origin, this results in an energy shift for 
S-wave states. 

In the same manner as for the chromomagnetic oper- 
ator, we can tune the coefficient C2 of the Darwin term 
so as to correct for the difference between the loop cor- 
rections of QCD and NRQCD. Previously, these correc- 
tions were unknown, and the coefficient C2 has been set 
equal to one in most non-perturbative simulations so far. 
While non-perturbative studies have shown a compara- 
tively mild dependence of energy levels on varying the 
value of C2 48J , a precise determination of the radiative 
corrections to c-i is a logical continuation of our improve- 
ment programme. 



A. Continuum QCD Calculation 

In continuum QCD, we need to consider the q 2 - 
dependence of the effective action, 

r« CD = *fi( 9 2 U* + + . . . , (36) 

2m 

where for the chromomagnetic case we had Fi(0) = 
Z\p(= Z 2 l ) and ^2(0) = Z a . Upon performing the non- 
relativistic reduction and isolating the terms containing 
the time component Aq of the gauge field, we find 



T® CD =F 1 (q 2 )^ \ 9 A 



_9_ 



2 q 2 A 



4m 2 



q 2 A 



(37) 



Ma 


1.9 2.65 3.4 


5?(i) , 7(1) , 7(1) 

Zj<j -+- Z> 2 T 
rytad , ( 1 ) | ryt ad , ( 1 ) 


-4.33(8) -4.16(7) -4.16(5) 

3.93 3.63 3.45 
0.78(8) 0.76(7) 0.83(5) 



TABLE III. Renormalization parameters for Chromomag- 
netic term, full spin v 6 NRQCD n = 4 



which, after renormalization, gives the Darwin term as 



[l-8m 2 F{(0) + 2F 2 (0)] ^ 



.9 2 A 

—2Q A 



■in i 



1pR + -- 



(38) 



Note that F 2 contributes through the non-relativistic 
reduction, whereas F\ only contributes through its ex- 
pansion around q 2 = 0, since the wavefunction renormal- 
ization cancels any contribution from Fi(0). Again, we 
do not have to include mass renormalization since the 
renormalized mass m naturally appears when using the 
non-relativistic reduction between on-shell spinors. 

Fa(0) has already been computed for the chromomag- 
netic case. For the contribution from the abelian diagram 
of Figure [ij a), one finds 



1 

(in 



9tt 



log /i/m 



(39) 



and for the non-abelian diagram of Figure [TJb) we may 
take the derivative of the analytical calculation carried 
out for the chromomagnetic term to obtain 



a. 



/J, 2 TT 



7m 
4^ 



11 

2ir 



9 



log /j,/m 



(40) 



The total contribution to the continuum QCD Darwin 
term is then 



(41) 



1 



fJ, 2 7T 



7m 1 , 6 
1 + 

4/1 7T 7T 



9/T 



I \ogfi/m 



where as before, we use a gluon mass fi as the IR regula- 
tor. We note, however, that the case of the Darwin term 
is more subtle since there are power-law IR divergences 
which will have to match and cancel with corresponding 
IR divergences on the NRQCD side. 



B. NRQCD Calculation 

The Darwin term in the NRQCD effective action is 
given by 



pNRQCD = _ Z NR gtft^q. 

D yv 8M 2 ^ 



(42) 
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where again Z|5 R has a tree-level contribution, Z™ 1 
1 + a a + . . .. After renormalization we obtain, 

pNRQCD = _ C2Z NR (C)2z NR # ,A fe + 

(43) 

Requiring the coefficients of the Darwin term to be equal 
in QCD and NRQCD, we find the matching condition 

c 2 ZTzr(ZT) 2 = Z D , (44) 

which at tree level and one-loop order gives 



c (0) - 1 

Jl) _ 7 NR,(1) 
c 2 — Zj d - Zj d 



- Z 2 NR ' (1) - 2ZN R ' (1) • (45) 



As in the chromomagnetic case, the diagrammatic con- 
tributions, labelled as Z, must be supplemented by the 
corresponding mean-field corrections, Z tad . 

Again, the symmetries of the effective action in BFG 
restrict the form of the diagrammatic contributions to 



Working in the Breit frame (pi = —qi/2) for ease of im- 
plementation, we isolate the renormalization constant 

7 NR,(1) _ 

AM 2 d 2 d 

Y,dq2 ST M(0,0)-M— ST^ Ao (0,0), 



dpo 



(47) 



where the second term arises from the on-shell condition 
of the incoming and outgoing quarks. We note that the 



final result is of course independent of the choice of frame 
used. 

By taking advantage of the modular structure of the 
HiPPy and HPSRC packages, we are able to reuse the 
same code that we used for matching the chromomag- 
netic operator by merely changing the incoming gauge 
field Lorentz index to isolate the A component, and tak- 
ing the trace of the diagram to isolate the implicit Dirac 
unit matrix in front of the Darwin operator. The pole 
structure for each individual diagram remains as in the 
previous calculations. 

The presence of severe IR divergences in some of the 
diagrams makes it necessary to analytically identify and 
subtract those divergences, leaving only the IR finite 
piece Z to be computed numerically. For the sake of 
brevity, the superscript NR is suppressed in the remain- 
der of this section except where it is necessary to avoid 
confusion. 

For the diagram in Figure [l|a) we obtain 



J D 



' 4 




9n ~ 


(err). 



log ^a, (48) 



where the IR divergence in round brackets arises from 
the inclusion of the tree-level Darwin term in the action. 
We therefore would expect this divergence to be cancelled 
by an IR divergence in the wavefunction renormalization 
times the tree-level operator. Since the IR divergences 

are calculated analytically we can use a constrained fit to 

?(«),(!) 



find the constant /ia-independent contribution, Z D 
For the diagram in Figure |T|(b) we have 



Z 



D 



D(Ma) log fia ■ 



7M 
4/7 



M 2 



(49) 



where D(Ma) is an action-specific coefficient containing 
lattice artifacts. The leading continuum power-law di- 
vergences are removed by the subtraction function 



I sub ( M ) = (4tt)(8M 2 ) 



d 4 k 

7M 21 

4fl 47T 



7 

8JU 



(3M 2 ) 



(k 2 + 
log /id 



H 2 ) 2 {ik + k 2 /2M) 
M 2 



{k 2 +n 2 ) 



2\4 



(50) 



where P sub (fia) is a polynomial in \xa with P sub (0) = 0. 
This subtraction function is chosen to cancel the leading 
M 2 / fj, 2 and M/fi IR divergences pointwise and also con- 
tains a con tinuum-like IR log divergence which is shown 
in Eq. |5o|. We calculate Z sub -W from a fit to I sub {^), 



using values of the cutoff in the range 10 < fx < 10 
and fitting the results to a polynomial in \ia. In Figure [4] 
we show a typical fit for the subtraction function giving 
the required renormalization constant Z sub . From the 
figure we note that the magnitude of the error increases 



as fi decreases and the cost of computation correspond- 
ingly increases. 

In addition to containing continuum IR divergences, 
the graph in FigureQb) contains additional artifact loga- 
rithmic IR divergences in NRQCD which must ultimately 
cancel against similar contributions from the graphs in 
Figure [ljc-f). We calculate these diagrams, namely the 
swordfish, algae, and ankh diagrams, and add their con- 
tribution to that of the subtracted diagram in Figurejljb) 
to compute the overall logarithmic IR divergence which 
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-3 -2 

log (//a) 



we compute, as a function of fia, the combination 

A(i) 



Z" 



jsub _ ^sub,{l) _j_ />(!) 



Z 6,(l) + Z c-/,(1) + ^_ log ^ a + p6-/ (/ia) (53) 

where P b ~'(/j,a) is a polynomial in /ia with P &_ ^(0) = 0, 
and the fit to the parameterization with the coefficient 
of the logarithmic IR divergence constrained to its ana- 
lytic value gives a much better estimate for the required 
renormalization constants. 

The presence of the tree-level Darwin term in the 
NRQCD action means we have also to include the con- 
tribution from the wavefunction renormalization 



,6,(1) 



Ma 



1.9. We obtain 



FIG. 4. Sample fit for Z£ 
-0.91(1) - 1.18(5)^0 with x 2 /d-o.f. = 0.21. 



should not (and indeed does not) contain a lattice arte- 
fact IR logarithm. 

For the swordfish, algae and ankh diagrams we write 



Z c D f ' {1) = Z c D J(1> +D'(Ma)\o gt Mi, (51) 



where D'(Ma) is the action-specific coefficient of the lat- 
tice artefact logarithmic IR divergence. We find that 



D(Ma) + D'(Ma) 



6 
- 



(52) 



which is independent of lattice artifacts, as expected. As 
before, the term in round brackets arises from the inclu- 
sion of the tree-level Darwin term in the NRQCD action. 
The graphs in Figure [TJb-f ) are combined with the sub- 
traction function in Eq. ( |50| and the integral computed 
numerically. Combining the diagrams in Figure [ljb-f) 



Hi) 



3tt 



log fia. 



(54) 



As expected, the IR logarithm in Z-i does indeed cancel 
the sum of the logarithmic IR divergences displayed in 
round brackets in Eqs. (48 1 and (52). 



The mass renormalization is finite and so can be simply 
added to the result. 

As for the chromomagnetic term, we must also include 
tadpole contributions, which in this case come from three 
sources: the improvement of the field strength tensor, the 
improvement of V^, and from cross-multiplication with 
tadpole-improved links in the other parts of the NRQCD 
action. Altogether, this gives 



SZ% d 



16 

y 



4M 3 32Af 2 



13 

AM 



,(2) 



(55) 



which is the same as the tadpole correction to the chro- 
momagnetic term plus 1 from the correction of \A ^ . The 
mass renormalization tadpole is the same as used in the 
chromomagnetic calculation. 

All lattice artefact IR divergences cancel internally 
within the NRQCD calculation and the continuum IR 
log divergences match between the continuum and lattice 
NRQCD calculations. Using the result for the continuum 
contribution in Eq. (41 1 and summing all lattice NRQCD 
contributions, we find the one-loop improvement coeffi- 
cient for the NRQCD Darwin term to be 



Jp = _I _ ^.(D _ - Z a >« - Z^ dA1) - Z^ - 2Z$ - 2Z% d 'W + 50 log Ma , (56) 



where again the logarithmic divergences in NRQCD and C. Results 

QCD have combined to give a term proportional to 

l °S Ma - l n Tables [TV] and [V] we give results for the full 



NRQCD action and the full spin v 6 NRQCD action re- 
spectively, both with stability parameter n = 4. For 
brevity, we use the notation 



(57) 
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Again, we find that the tadpole correction cancels the 
major part of the diagrammatic contribution to Cq , in- 
dicating improved convergence of the resulting pertur- 
bative series. This effect is expected and reinforces the 
efficacy of including tadpole improvement. 

Numerical simulations 48J show that using a value of 
C2 = 1.25 shifts the energy of the T(15) state by less than 
one percent compared with the tree-level value c 2 = 1, 
indicating that the radiative corrections are well under 
control. 



Ma 






(!) 

c 2 av C2 


1.9 


-5.95(8) 


5.18 


1.33(8) 0.22 1.29(2) 


2.65 


-3.71(10) 


4.88 


0.08(10) 0.25 1.02(3) 


3.4 


-1.73(12) 


4.69 


-1.20(12) 0.27 0.68(4) 



TABLE IV. Renormalization parameters for Darwin term, full 
NRQCD n = 4 



Ma 


Za-f,(\) 


rrtad,(l) 


(i) 

c 2 ay C2 


1.9 


-5.95(8) 


5.18 


1.33(8) 0.22 1.29(2) 


2.65 


-3.71(10) 


4.88 


0.08(10) 0.25 1.02(3) 


3.4 


-1.73(12) 


4.69 


-1.20(12) 0.27 0.68(4) 



TABLE V. Renormalization parameters for Darwin term, full 
spin v° NRQCD n = 4 



VI. FOUR-FERMION INTERACTIONS 

At the one-loop level, calculations in NRQCD must 
also take account of four-fermion interactions in the 
NRQCD action that are needed to match QQ — > QQ 
scattering processes between QCD and NRQCD. While 
so far we have emphasized the role of the effective action 
in BFG, the matching condition must ultimately equate 
physical, on-shell matrix element, which in general will 
include contributions from 1PR diagrams, and so it is 
not necessarily possible to match just the effective ac- 
tions term by term. 



A. Formalism 

The radiatively generated four-fermion interactions 
that need to be added to the NRQCD action can be 
written as four-fermion contact operators in a (covari- 
ant) derivative expansion. Here we consider only the 
lowest-order terms in this expansion, namely those with- 
out derivatives. One obvious set of operators which bear 
a close relationship to the QQ — > QQ diagrams calcu- 



lated, is given by 

gNRQCD _ 



«8 



M 2 



(X c*x)(-0o"0) • 
(58) 



We will refer to the operators with coefficients a\ and b% 
as singlet-exchange operators. The operators with coeffi- 
cients as and &8, which we will refer to as octet-exchange 
operators, give corrections to processes involving single 
gluon exchange at tree level. This relation will be useful 
later when discussing corrections to the QCD Coulomb 
force. The spin-independent and spin-dependent oper- 
ators have coefficients ai and bi, i = 1,8, respectively. 
Only b\ and bg will contribute to the hyperhne splitting, 
as can be seen by choosing the more conventional set of 
operators given by [H [5] 

gNRQCD = di A^ x * ){x T^ + da^WW**) 

+d 3 ^(^T aX *)(x T T a ^)+d 4 ^(^aT aX *)(x T ^) . 

(59) 

Here, we choose the X field for the antiquark to transform 
according to the conjugate 1/2 representation of spin and 
the 3 representation of colour generated, respectively, by 
— a* and —T*. Such s-channel operators make explicit 
their effects on different meson states; this will be useful 
in later discussions of the hyperfine splitting. 

These two sets of operators in Eqs. (58), (59) are re- 



lated by Fierz transformations. Considering the action of 
the operators on colour octet (-^T b ) and singlet states 

(-L/g), and spin-0 {^h) and spin-1 (-Ter 3 ) states, we 
find the relevant Ficrz transformations to be 



di = An 

d 2 = 47T 
d 3 = 47T 
d 4 = 47T 



46 8 + gfls + 36i + ax 

4 , 4 
3 1 

7^8 - ^a 8 + 3&i + ai 
2 6 



\b 8 - \a 8 -h+ai 
2 b 



(60) 



Note that the factor of 47r is needed to re-express the 
contact operators in terms of a 2 s rather than a s g 2 . 
From Eq. ( 59 ) we see that the energy of colour singlet 



QQ mesons receives a contributions proportional to d\ 
in the spin singlet case and to d 2 for the spin triplet 
case. The hyperfine splitting is therefore proportional to 
(d\ — d 2 ), cf. subsection VII|VII A for further discussion. 
From Eq. ( 60 ) we then find 



16tt 2tt 

~T b8 + T ' 



(61) 
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(a) (b) 




(c) (d) 

FIG. 5. Ladder Feynman diagrams. Referred to as (a) box, 
(b) crossbox, (c) triangle, and (d) diamond diagrams. Note 
that there are two diagrams with the triangle topology. 



pairs. Figure [6] shows the relevant continuum QCD dia- 
gram that is absent in NRQCD. Labelle et al. [8] give the 
amplitude for this process for QED and it is straightfor- 
ward to obtain the corresponding QCD result by includ- 
ing the correct colour factors to give [25] 

d 2 -"' ann = - ?|» (2 - 2 In 2) . (62) 



and so the coefficients a\ and as of the spin-independent 
operators have cancelled out. In what follows we there- 
fore calculate the coefficients 6g and bi of the spin- 



dependent operators in Eq. ( 58 1 



Since there is no <x • A coupling in NRQCD, it turns 
out that 1PR diagrams cannot contribute to the spin- 
dependent non-derivative operators and hence the oper- 
ator combination which gives rise to the hyperfine split- 
ting can be matched using 1PI diagrams only. In or- 
der to complete the matching of all four-fermion oper- 
ators, including the spin-independent contributions, we 
need to also calculate the 1PR contributions from insert- 
ing the vacuum polarization onto a Coulomb gluon ex- 
change line, since this will generate four-fermion contact 
operators to all orders in the derivatives. It also turns 
out that the IR divergences in the spin-independent case 
are more severe than for the spin-dependent case due 
in part to these contributions from Coulomb exchange, 
which also require a careful treatment of lattice artifacts 
multiplying logarithmic IR divergences. We shall com- 
ment further on this observation below and will present 
the general matching analysis for all non-derivative four- 
fermion operators in a future paper. 

As far as the spin-dependent interactions are con- 
cerned, the leading contributions to the four-fermion in- 
teractions arise from the 1PI ladder diagrams shown in 
Figure [5j The box (a) and crossbox (b) diagrams con- 
tribute in both continuum QCD and NRQCD, whereas 
the triangle (c) and diamond (d) diagrams are specific 
to NRQCD; again, let us emphasize that these NRQCD- 
specific diagrams consist not only of lattice artifacts, but 
also contain continuum contributions from vertices aris- 
ing from the non-relativistic form of the NRQCD action. 
Much as in the previous sections, by considering the 1PI 
diagrams in both continuum QCD and lattice NRQCD, 
we can then determine the values of the coefficients of the 
spin-dependent four-fermion operators in the NRQCD 
action needed to account for radiative improvement cor- 
rections. 

In addition to the corrections due to QQ — > QQ scat- 
tering calculated above, a further correction to d\ is re- 
quired to account for the fact that the NRQCD formal- 
ism does not allow for the creation or annihilation of QQ 




FIG. 6. Two-gluon annihilation diagram 



B. Calculation 



The absence of a QQ vertex in NRQCD and the fact 
that quarks and antiquarks have identical masses by the 
CPT theorem means that the simplest implementation of 
antiquarks using the HP SRC package is to treat them as 
just a different species of quark. In fact, we can simply 
relate QQ-scattering diagrams 



2 

+ Z^ll-^(i)lia(iT a iT b ) a $i)$) (^ia(iT a iT b )/3^) 

2 



(63) 



to the QQ one by changing the representation of one of 
the quarks appropriately. For NRQED, the antiparticle 
vertex can be obtained by replacing e — > — e; the NRQCD 
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analogue is taking (a, T a ) 



-a*, -Tj) to obtain . 



- z ^-l^p(xi^*(T b T a )l s x8)(^l(T(T a T b )p 7 4> 7 ) 



- z Zm-a^r(xi^*(T a T b )l 5 xs)(^l<y{T a T b )p^) , 

(64) 

where the transpose of the individual T a elements has 
been pulled outside and the antiquark fields have been 
relabelled to \ to agree with convention. We have kept 
all indices explicit to highlight how the creation and an- 
nihilation operators are directly replaced. 

Dividing out the colour factors and isolating the cor- 
rect spin structure on both fermion lines passing through 
the diagram, we can calculate the one-loop coefficients, 
Zfym etc., directly in the HPSRC code. The factors of 
i in Eq. ( 63 Icome from the convention that we work with 
anti-hermitian generators when evaluating the colour fac- 
tors of diagrams. This is convenient since then the com- 
mutation relations have real structure constants. Note 
that by construction the box diagram (a) contributes 
only to the symmetric colour combination, whereas the 
crossbox diagram (b) contributes only to the antisym- 
metric colour combination. The triangle (c) and di- 
amond (d) diagrams contribute to both colour combi- 
nations since they contain two-gluon vertices carrying 
colour factors {T a ,T b }. 

Although the momentum exchange acts as an IR regu- 
lator in these diagrams, we set all three- momenta to zero 
with the external quarks on-shcll and use a gluon mass fi 
as the IR regulator using a Stiickelberg mass term. This 
is possible as discussed above, and helps the evaluation 
by significantly simplifying the pole structure of the in- 
tegrand: all diagrams have a simple pinch at fco = ±/i, 
and the calculation can be performed without needing to 
implement contour shifts in the complex energy plane. 
As in the case of the BFG vertices, the absence of UV 



divergences allows us to employ different UV regulators 
in the QCD and NRQCD calculation. 

Radiative improvement is then achieved by adding 
four-fermion operators to the lattice NRQCD action 
with coefficients tuned such that one-loop calculations 
of QQ — > QQ scattering in NRQCD give the correct con- 
tinuum QCD result. For example, we add terms such 
as 

2 

r = -Zilln^^ixUTbTa^sXs)^^^)^^) , (65) 

where ziy m — Z^m^ — zfy^n^ is the difference be- 
tween the diagrammatic coefficients calculated in contin- 
uum QCD and NRQCD, respectively. The radiative im- 
provement can then be implemented with the coefficients 
given by 



as = 


+ 1)411 + (- 


-3 + 


-)z<n i 

<r^J asym) ) 


(66) 


ai = 


_ 2 (^(i) + zW I 

q I sym 1 asymS 5 






(67) 


h = 




(-3 


4- 5 )7 {1) 

' ^ ' asym- 


-.1,(68) 


h = 


_ 2 r 7 (i) , 7 m 

g X syrn — a ' asym- 


-J- 




(69) 



The analytical calculation of the continuum QCD dia- 
grams is relatively simple with our choice of IR regulator. 
The results are 



7 



Z QCD,(1) 

sym " 12tt ' 4tt 



l / i \ 5M M 2 M 3 
-IogOi/M) + — - 7 ^ ~*~ 



^g MM) ^ , ^ 



16/i 2tth 2 ' 2/i 3 
3M M 2 



■QCD,(l) 
sym— a 

■QCD,{1) 
sym — <j 



— + _log( M /A/) + - 
1 



-ilT 



log (fi/M) 



(70) 



(71) 



The NRQCD diagrams are computed numerically us- 
ing the HPSRC package. Since higher order IR diver- 
gences dominate numerically, suitable subtraction func- 
tions must be chosen to pointwise cancel the divergent 
integrands: 



r ub (u) — 

sym VrV 


-/ 


r ub (u) — 

asym VrV 


-/ 


I sub (a) — 

1 sym—crKH 1 ) 


-/ 


I sub (a) — 

asym— a VrV 


1 


2M 2 



(1 + fc 2 /4M 2 ) _ ~ h _ I 5M M 2 _ M 3 

\k 2 + M 2 ) 2 (fc 2 + fc 4 /4Af 2 ) ~ Z «f»W 7T i0g [ ^ a) 16/i + 2nn 2 2y? 



(72) 

1 ~ 15 3M M 2 

d4k (k 2 + v 2 ) 2 (*h + k 2 /2M) 2 = z -y> M + 8^ log ( ^ a) + ^7 - 2^ ' (73) 

k 2 12M 2 ~ 1 M 

dAk wrwm^i^) = Zlt - M ~ * log M " ■ (74) 

d ' k jWTW? = ^ym-M + ^ log M • (75) 



i 

These subtraction functions are continuum-like in that they contain no lattice artifact IR divergences. Evaluat- 
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-0.08 
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-6 




-5 



-4 -3 
log (/id) 



-2 



-1 



Ma 


1.95 


2.8 


4.0 


6i 


-0.0080(3) 


-0.0319(4) 


-0.0608(5) 


b 8 


0.0274(8) 


-0.0436(15) 


-0.1451(26) 



TABLE VI. Renormalization parameters for spin-dependent 
four-fermion operators, full NRQCD n = 2 



Ma 


1.9 


2.65 


3.4 


bi 


-0.0043(1) 


-0.0266(6) 


-0.0459(2) 


b 8 


0.0378(5) 


-0.0266(1) 


-0.0909(8) 



TABLE VII. Renormalization parameters for spin-dependent 
four-fermion operators, full NRQCD n = 4 



FIG. 7. Sample fit for Z™^_ CT (/i) for Ma = 1.95. We obtain 
-0.0159(17) - 0.145(15)^a + 0.08(3)(mci) 2 with ^jd.o.f. = 
0.09. 
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FIG. 8. Sample fit for Z s a t ym _ a {\i) for Ma = 1.95. We obtain 
-0.186(3) + 0.3(l)/ia - 0.2(3)(>a) 2 with xV d -°-/- = 0.08. 



ing I suh numerically and fitting (with the coefficients of 
the divergences constrained to agree with the analytical 
results, and a polynomial in \ia added) allows us to deter- 
mine Z sub . As the integrands are easily evaluated, a large 
number of points can be sampled. The high-order diver- 
gences cause the integrand to become extremely large 
around k ~ 0, which can mean that for small enough 
k even double-precision variables can overflow, thus re- 
turning NaN for the value of the integrand. To avoid this, 
we add a statement that will set the integrand to zero 
in a small neighbourhood of the origin. We have to use 
a sufficient number of VEGAS sampling points in order 
to resolve this cut around the origin exactly. In Figures 
([7])-([8]) we plot sample fits of these subtraction functions. 

For each colour ordering, the sum of all diagrams is 
then calculated together with the appropriate subtrac- 
tion function. In this way all IR divergences except, in 
some cases, the logarithmic IR divergence are cancelled, 
regardless of the lattice artifacts, allowing a more con- 
strained fit. For the spin-dependent contribution we cal- 



culate 



^ syra — a ' * sym- 



7 NR,(1) 



\ -'-as 



ct (m) 
00- 



*-* sym- 
zsub 



, m _ ~NR,(1) 
T \yj) ^sym — a 

/n\ _ 5fNR,(l) 
ry^/ ^asym — c 



(76) 
(77) 



We evaluate the diagrams for several gluon masses in 
the range 10 -3 < /i 2 < 1CP 2 , which is chosen such that 
any lattice artifacts (appearing as polynomials in fia) are 
negligible, yet at the same time large enough such that 
we do not run foul of any hard- wired numerical tolerances 
in the HPSRC package. 



Results 



We then have 

syrn — a 
asyrn — a 



?NR,(1) 
J sym—o 

;nr,(i) 



l 

127T 
1 



-in 



log M a , 



Using these equations, the numerical results for the ra- 
diative improvement coefficients b\ and 6g for the spin- 
dependent operators are calculated from Eq. (68) and 
( |69[ ). Results for various NRQCD actions are given in 
Tables (VI I to (VIII). Note that the coefficients have a 
sizeable dependence on Ma. For the full spin v 6 NRQCD 
action, this dependence is less pronounced. 

As noted at the beginning of this section, we will re- 
port on the calculation of the radiative improvement co- 
efficients a\ and ag for the spin-independent operators in 
Eq. ( 58 ) in a future paper. 



Ma 


1.9 


2.65 


3.4 


bi 


-0.0341(2) 


-0.04778(2) 


-0.0622(2) 


bs 


-0.0899(9) 


-0.1243(7) 


-0.1690(9) 



TABLE VIII. Renormalization parameters for spin-dependent 
four-fermion operators, full spin v 6 NRQCD n — 4 
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A. Hyperfine Splitting 



FIG. 9. The tree-level diagram contributing to heavy-heavy 
meson hyperfine splitting: a spatial gluon is exchanged be- 
tween two vertices involving the chromomagnetic operator. 



VII. APPLICATION TO T AND r) b SPECTRUM 



Having derived the radiative corrections to the a ■ B, 
Darwin, and spin-dependent four-fermion operators in 
the NRQCD action, we proceed to analyse the effects of 
these corrections on the spectrum of mesons containing 
b-quarks. Depending on whether the operators involved 
are spin-dependent or not, we can distinguish between 
changes to the hyperfine splitting and overall shifts of 
the ground state mass. While obviously no substitute 
for simulations including the radiatively corrected coeffi- 
cients in the NRQCD action, estimates for these effects 
are important because they give a clear indication about 
the expected magnitude of the corrections to the spec- 
trum of S-wave bb states. In what follows the quark mass, 
M, is identified with the mass, Mb, of the &-quark and 
defined on the lattice as a~ 1 {aM b ). 



The leading contribution to the hyperfine splitting can 
be estimated from a perturbative picture involving the 
exchange of a single gluon between two vertices involving 
the chromomagnetic operator (cf. Figure [9]); this process 
shifts the energy of a colour-singlet meson state \M) = 
rf + bb + gg) \ S) by 



l 



AE: 



(M\(axtqT a )(a* xiqT^)\M}/q 2 



4M 2 

~{S\(v x iq)(a* x lq )\S) l -Tr{T a Ta) / 'q 2 



-241 

9M 2 



(S\a.a*\S) . 



(78) 



For \r) b ) = -±\ 1 1+ I 1) (S=Q) and |T >= 1 1 1) (S=l), 
we have (r] b \a.a* \r] b ) - (T\a.a*\T) = 3 - (-1) = 4, and 
hence the hyperfine splitting is approximately 



M, 



9M 2 



1-0(0)1- 



(79) 



where ip(0) is the meson wavefunction at the origin, as- 
sumed to be the same for the r)h and T states. We see 
that the rjb lies below the T, as expected. 

Nonpcrturbative results indicate that the hyperfine 
splitting in heavy quarkonium is indeed approximately 
proportional to c 2 [S]. Given that the origin of these 
nonperturbative results can be understood from a tree- 
level estimate, we can proceed to examine the effect of 
the four-fermion operators in the same manner. This is 
straightforward when working in the basis of operators 
given by Eq. ( 59 1 , where that we get 



M„. - M r 



9(dr - d 2 ) 4 
2 3M 2 



1-0(0)1- 



(80) 



for the four-fermion contribution to the hyperfine split- 
ting, which can be re-expressed in term of the coefficients 
Oi, bi using Eq. (60). 



The one-loop correction to the hyperfine splitting is 
then estimated to be the tree-level value multiplied by 
the one-loop correction factor 



1 + a v (q*)[2c 



9 fl6 

-8 U 68 



4&i 



8tt 



(2 -2 In 2) 



(81) 



In Table IX 



we give our estimates for the corrections 
that need to be applied to the hyperfine splitting as mea- 
sured in full v 4 NRQCD n = 2 without perturbative im- 
provements. Applying our estimated correction post hoc 
to the data presented by the HPQCD collaboration in [1], 
we find that the corrections to the chromomagnetic oper- 
ator and the inclusion of the spin-dependent four-fermion 



terms work to increase the hyperfine splitting, pushing it 
closer to the experimental value of 69.3(2.8)MeV [lOj . 
We note that the corrections from the four fermion op- 
erators further reduce the lattice spacing dependence, as 
can be seen from the results plotted in Figure 10 It is 



important to note that for the 0{v A ) NRQCD action the 
errors on the corrected results include an estimate of the 
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M p ay 


C4 correction box correction total 


old Ms 


new hfs 


1.95 7.09 0.216 
2.8 6.76 0.249 
4.0 6.458 0.293 


+31.4(3)% -4.4(1)% +27.0(3)% 
+39.8(3)% +8.3(2)% +48.1(4)% 
+49.3(3)% +31.3(3)% +80.9(4)% 


56(2) 
50(2) 
41(2) 


71(3)(5)(6) 
74(3) (6) (5) 
74(3) (7) (4) 



TABLE IX. Estimates of the corrections to the bottomonium hyperfine splitting results of [TJ arising from the radiative 
improvement of the n — 2 full v 4 NRQCD action. The errors given in the last column are statistical, 0(a 2 ), and relativistic 
corrections, in that order. 



o 


Uncorrected 


n 


Corrected 


o 


PDG 



a" (fm) 



FIG. 10. Comparison of the corrected and uncorrected bot- 
tomonium hyperfine splitting results for n = 2 full v 2 NRQCD 
[I]. Note that for the corrected results the total error includ- 
ing estimates of the effects of 0(ct 2 ) contributions and 0(v 6 ) 
terms omitted in the simulation is displayed, whereas the un- 
corrected results are shown with statistical errors only, since 
their 0(a 3 ) errors would be too large to show on this scale. 
The PDG data point is taken from |10j . 
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FIG. 11. Comparison of the corrected and uncorrected bot- 
tomonium hyperfine splitting results for n = 4 full v 4 NRQCD 
[48] . Note that for the corrected results the total error includ- 
ing estimates of the effects of 0(a 2 ) contributions and 0(v 6 ) 
terms omitted in the simulation is displayed, whereas the un- 
corrected results are shown with statistical errors only, since 
their 0(a s ) errors would be too large to show on this scale. 
The PDG data point is taken from 10 • 



Ma ay(q*) 


hfs (MeV) 
C4 = 1 improved C4 


Correction 
4-fermion 


hfs (MeV) 
corrected 


1.9 0.225 
2.65 0.253 
3.4 0.275 


56.1(1) 72.1(1) 
50.5(1) 69.8(1) 
45.6(1) 65.6(1) 


-1.7(1) 
+4.5(1) 
+10.3(1) 


70.4(1)(28)(56) 
74.3(1)(32)(50) 
75.9(1) (34) (46) 



TABLE X. The corrections to the bottomonium hyperfine 
splitting arising from the radiative improvement of the n = 4 
full v 4 NRQCD action, as found in [48] . Only the four-fermion 
contributions are post hoc estimates. The errors given in the 
last column the errors are statistical, 0(a 2 ), and relativistic 
corrections, in that order. 



effect of the omitted 0(v 6 ) term shown in Eq. (J8|, and 
that these errors are preliminary in the sense that future 
simulations will include these terms explicitly, using the 
results of this paper, and removing the need for a post 
hoc correction. 

In Table [X] and Figure [IT] we present the corrections 
to the full ir NRQCD n = 4 hyperfine splitting as mea- 
sured in [48], where the perturbatively corrected coeffi- 
cient C4 has been included in the simulation, and only the 
four-fermion operator corrections need to be applied by 



M olv 


C4 correction box correction 


1.9 0.225 
2.65 0.253 
3.4 0.275 


+37(4)% +17.3(1)% 
+40(4)% +26.2(1)% 
+26(4)% +37.6(1)% 



TABLE XI. Corrections to the bottomonium hyperfine split- 
ting arising from the radiative improvement of the n = full 
NRQCD action including spin-dependent terms at order v 6 . 



hand. As in n — 2 case, the corrected chromomagnetic 
operator acts to increase the hyperfine splitting, and the 
four-fermion corrections reduce the lattice spacing depen- 
dence (which however is less severe from the outset when 
comparing to the n — 2 case). 



Finally, in Table XI we estimate the relative corrections 
from the corrected coefficients to full n — 4 NRQCD re- 
sults when including the spin-dependent v 6 terms in the 
action. While the corrections arising from the improve- 
ment of the chromomagnetic term are similar to the case 
of full n — 4 v 4 NRQCD, the corrections coming from the 
four-fermion operators are significantly larger. While one 
might naively expect that the effect of radiative correc- 
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tie a 2 dependence, and even for coarse lattices are in 
very good agreement with experimental data; this gives 
a strong check on the correctness and usefulness of the ra- 
diative corrections and lends credibility to the prediction 
(rather than postdiction) of Mb- - M Bc = 54(3) MeV 
|51j incorporating the effects of radiatively improving C4. 



B. Mass shift 



FIG. 12. The tree-level diagram contributing to the hyperfine 
splitting in heavy-light mesons: a spatial gluon is exchanged 
between two vertices, the heavy-quark one of which involves 
the chromomagnetic operator. 



tions ought to decrease as higher-order terms are added 
to the action, it has been shown in [5] that the inclusion of 
0(v 6 ) terms (albeit with a different gauge action) leads 
to a decrease of the (tree-level) hyperfine splitting, ft 
appears that the larger corrections from the four-fermion 
operators would compensate for this effect. 

The radiative corrections to the chromomagnetic oper- 
ator (but not the four-fermion operators) will also affect 
the hyperfine splitting in heavy-light mesons through the 
leading perturbative contribution shown in Figure 12 
giving a hyperfine splitting approximately proportional 
to C4. Besides the absence of the four-fermion terms, 



heavy-light systems also benefit from much smaller v 
corrections, which makes them a particularly suitable 
test case for assessing the efficiency of including the ra- 
diative corrections. In [51], the NRQCD action with ra- 
diatively improved coefficients was used for the first time 
for the heavy-light £?-meson system. The results for the 
hyperfine splittings for the B and B s mesons show lit- 



The leading spin-independent perturbative correction 
to the energy of a meson state is given by the single- 
gluon exchange involving the Darwin term at one of the 
vertices, as shown in Figure |13| This gives the energy 
shift from the corrected C2 coefficient as 



2c 2 g 2 
8M 2 
-C29 2 
MI 2 



(M\(q 2 T a )(Tj)\M)/q 2 

\m\ 2 - 



AE* ^4IL v {U\u?T„HT. J , ) U)/« 2 (82) 

(83) 

The corresponding contribution from the four-fermion 
operators gives 



AE 



9rfi 4 
~2~3M 2 



WW 



(84) 



While the hyperfine splitting between the rjb and T 
states will lead to higher-order corrections to their wave- 
functions at the origin, to leading order these can be 
taken to be identical, here denoted -0(0) , giving the same 
shift in mass for both states. 

Again, there is no real substitute for including the four- 
fermion operators in numerical simulations; however, we 
can attempt to estimate the size of their overall effect by 
determining an effective value for c 2 using the tree-level 
approximations given above. 



Combining Eqs. (83) and (84), we find 



AE 



3M 2 



1 + a s (c ( 2 > - I2b s - 4a 8 - %i 



3ai + -(2-21n2)) |V>(0)| 5 

7T 



(85) 



giving an effective value of = 1 + a s (c^ — 126 8 — 
Aa% — 9&i — 3ai + - (2 — 2 In 2)). From this expression it 
is clear that we cannot consider the correction due to the 
Darwin operator in isolation, but must include the effects 
of the spin-independent four-fermion operators which we 
have not computed in this paper and whose calculation 
will be presented in a forthcoming paper. 



VIII. CONCLUSIONS 

In this paper, we have applied the BF method to lat- 
tice NRQCD and have computed the one-loop radiative 



correction to the coefficient, C4, of the er • B and the one- 
loop radiative contribution to the coefficients, d\ and d 2 
of the four-fermion contact operators that affect the hy- 
perfine structure of heavy quark mesons. The gauge in- 
dependence of our calculation was explicitly checked by 
carrying out both relativistic and non-relativistic calcu- 
lations in the lattice theory. This is possible because in 
BFG all calculations are UV finite. 

Our results are summarized in Tables |IX| |Xj and |XI| 
and in Eqns. (35), (81) and (85). In particular, in eqn. 



( 35 1 there is a negative correction to C4 due the the effect 



of the continuum logarithmic IR divergence. However, it 
turns out that the constant terms more than cancel this 
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effect and the correction to C4 is positive. 

Whilst no substitute for including these corrections in 
a simulation, our estimate for the correction to the T — rjb 
hyperfine splitting measured by Gray et al. [T], as shown 
in table |IX| indicates that the effect of the corrections is 
to reduce the lattice spacing dependence to within the 
remaining errors. The resulting estimate for the hyper- 
fine splitting of 72(3) (5) (3) McV is then in good agree- 
ment with the experimental value of 69.3(2.8) MeV [ID] . 
Subsequent simulations [4*8l I5T] have confirmed these ex- 
pectations. It will be interesting to see how the inclusion 
of our radiative corrections in simulations utilizing the 
0(v 6 ) NRQCD action will compare with the results of 
[5], where a reduced hyperfine splitting was found when 
using the tree- level 0(v 6 ) NRQCD action; from our re- 
sults, we expect that the inclusion of the four-fermion 
operators will compensate for this reduction. 

The elimination of 0(a s a 2 ) errors, the much reduced 
dependence of observables on a 2 and the agreement with 
experiment gives us confidence that the improvement 
strategy for constructing the NRQCD effective action is 
robust. 
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FIG. 13. The tree-level diagram contributing to the ground 
state mass shift: a temporal gluon is exchanged between two 
vertices, one of which involves the Darwin operator. 
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